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ABSTRACT 

We compare the structural properties of galaxies formed in cosmological simu- 
lations using the smoothed particle hydrodynamics (SPH) code GADGET with those 
using the moving- mesh code AREPO. Both codes employ identical gravity solvers and 
the same sub-resolution physics but use very different methods to track the hydrody- 
namic evolution of gas. This permits us to isolate the effects of the hydro solver on 
the formation and evolution of galactic gas disks in GADGET and AREPO haloes with 
comparable numerical resolution. In a matching sample of GADGET and AREPO haloes 
we fit simulated gas disks with exponential profiles. We find that the cold gas disks 
formed using the moving mesh approach have systematically larger disk scale lengths 
and higher specific angular momenta than their GADGET counterparts across a wide 
range in halo masses. For low mass galaxies differences between the properties of the 
simulated galaxy disks are caused by an insufficient number of resolution elements 
which lead to the artificial angular momentum transfer in our SPH calculation. We 
however find that galactic disks formed in massive halos, resolved with > 10^ parti- 
cles/cells, are still systematically smaller in the GADGET run by a factor of ~ 2. The 
reason for this is twofold: i) the excessive heating of haloes close to the cooling radius 
due to spurious dissipation of the subsonic turbulence in GADGET reduces the supply 
of gas which can cool and settle onto the central disk; ii) the efficient delivery of low 
angular momentum gaseous blobs to the bottom of the potential well results in the 
centrally concentrated gas disks in GADGET simulation. While this large population 
of gaseous blobs in GADGET originates from the filaments which are pressure confined 
and fragment due to the SPH surface tension while infalling into hot halo atmospheres, 
it is essentially absent in the moving mesh calculation, clearly indicating numerical 
rather than physical origin of the blob material. 
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1 INTRODUCTION 

A primary goal of cosmological simulations is to self- 
consistently reproduce the variety of galaxy morphologies 
observed in the local Universe. While the formation of dark 
matter haloes via gravitational collapse has been simu- 
lated in great detail using N-body simulations (e.g.,'Springel| 
et al.||2005||Boylan-Kolchin et al.pOOQ; Fo salba et al. 20Q8| 
Teyssier et al.||2009| [Klypin et al.j|2011| ), modeling the evo- 
lution of the luminous components of galaxies has lagged 
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behind due to the intrinsic complexity of gas dynamics and 
star formation. Early efforts to incorporate baryonic pro- 
cesses into cosmological simulations accounted for gas cool- 
ing, but did not include star formation or related feedback 
effects. These studies found that cooling gas accreted into 
dark matter haloes would quickly lose angular momentum 
and fall to the centre of the potential ( Navarro Benz|1991 



Katz Gunn|| 1991 1 [Navarro White||1994t . The forming 



objects had disk-like morphologies, but with low specific 
angular momenta, and with most of the gas residing in a 
central spheroid rather than a rotationally supported disk, 
unlike most observed late- type galaxies. It was argued that 
the efficient angular momentum loss was largely a conse- 
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quence of the early collapse and formation of proto-galactic 
clouds which were able to efficiently transfer their angular 
momentum to the dark matter haloes by dynamical and hy- 
drodynamical friction during merger and accretion events. 

A large number of subsequent disk formation studies 
have attempted to diagnose and fix this so-called "angu- 
lar momentum catastrophe". Most proposed solutions are 
centered around preventing the gas from cooling and form- 
ing stars too efficiently at high redshift. In a simple exper- 
iment, Weil et al. (1998) showed that if gas cooling is pre- 



vented until 2; = 1, stellar disks could form with specific 
angular momenta consistent with observed spiral galaxies. 
The two widely advocated mechanisms to mitigate gas over- 
cooling are heating by an ultraviolet (UV) radiation field 
and feedback associated with star formation. The UV back- 
ground (Quinn et al. 1996; Navarr o Steinmetz|1997 Hoeft 
et al.|[2bQ6 ) has been shown to inhibit the accretion of cold 
gas by low mass haloes, but does not appear to provide a full 
solution to the angular momentum problem. Star formation 
with associated feedback has also been identified as a heating 
mechanism which might prevent early collapse (e.g.lThacker 



Couchman|2000[ [2001] [Mailer Dekel|2002| |Abadi et al. 
20031 IRobertson et al.||2004| [Qkamoto et al.||2005| |Scanna- ' 



pieco et al.||2008| ) or efficiently remove low angular momen- 
tum material (e.g. Governato et al.|2010 Gued es et al.|2011 ) 
which allows less centrally concentrated disks to form. Al- 
though strong feedback can improve galactic disk formation, 
it is not immediately clear that this is the only solution to 
the angular momentum problem or if other numerical arti- 
facts remain adversely affecting the formation of rotationally 
supported galaxies. 

There are some well studied issues with the standard 
density formulation of SPH - which is the most com- 
monly employed SPH formulation for cosmological simula- 
tion codes - that can cause spurious angular momentum 



transfer from gas disks. Okamoto et al. (2003) showed that 



gas disks embedded in a diffuse hot halo would systemati- 
cally loose angular momentum due to spurious hydrodynam- 
ical torques, an effect that is particularly severe at low reso- 
lution (see also Commercon et al.^2008j. However, t his prob- 
lem is resolution dependent, and Governato et al. (2004) il- 



lustrated this point by presenting disk galaxy formation sim- 
ulations in a ACDM context - without invoking strong feed- 
back - to show that angular momentum loss could be sub- 
stantially reduced by increasing the mass and spatial resolu- 



tion. Similarly, Kaufmann et al. (2007) used idealized inside- 



out disk formation simulations to show that while spuri- 
ous hydrodynamical angular momentum loss dominates at 
low particle resolutions, using > 10^ SPH particles in each 
simulated galaxy can make the unphysical hydrodynamical 
torques subdominant. 

Unfortunately, the very high resolution criteria speci- 
fied in Kaufmann et al. (2007) make the near-term feasibil- 



ity of carrying out full cosmological box simulations with 
standard SPH poor. As a result, many recent galaxy forma- 
tion studies have adopted the "zoom-in" technique, where a 
single galaxy can be simulated at a very high mass and spa- 
tial resolution (e.g. Guedes et al.||20lT Agertz et al.||2011 ). 
Some of these efforts have led to the formation of galax- 
ies that share many properties in common with our own 



ducing a Milky Way type disk galaxy - the properties of 
their simulated galaxies depend heavily on the choice of the 
star formation threshold, formation efficiently, feedback pa- 
rameters, and other poorly constrained star formation re- 
lated parameters. Since "zoom-in" simulations are limited 
in their scope to one halo at a time, it becomes difficult to 
judge if the same simulation parameters (e.g., star formation 
threshold, etc.) would validly reproduce the wide range of 
observed galaxy morphologies or observationally constrained 
quantities such as the global star formation rate at differ- 
ent redshifts. So, while "zoom-in" simulations are a very 
useful numerical tool to understand how individual galaxies 
form and evolve, it is necessary to perform large cosmological 
box simulations, where a wide variety of structures should 
naturally form which can then be compared directly to the 
wealth of observational galaxy data. Simulations of represen- 
tative samples of the Universe permit us to test more clearly 
the impact of poorly constrained simulation parameters on 
structure formation by exploring the evolution of a full en- 
semble of galaxies, rather than one individual object. Once 
it becomes feasible to produce a large ensemble of realistic 
galaxies in a cosmological context, we will have a powerful 
tool to address questions about the driving forces behind 
galaxy morphological evolution that would complement the 
efforts of "zoom-in" simulations. 

Another issue that exists in standard density SPH is the 
formation of dense g as "blobs" |K aufman n et al. || 2006 1 [van 
de Voort et al.| 2011| [Keres He rnquist 2 009| [van de Voort 



Schaye|2012 ) , that form via numerical thermal instability 



that occurs in the absence of thermal conductivity (Hobbs 
|et al.|2012). For example, |Kaufmann et al. (2006) presented 
simulations of inside-out disk formation and found a popula- 
tion of dense gas blobs efficiently accreted onto their central 
forming galaxy. These blobs - which are not found in adap- 



tive mesh refinement simulations { Joung et al. 2012| ) or more 
modern SPH algorithms where entropy mixing is included 
via thermal conductivity ( Hobbs et al.|2012 ) - deliver a sub- 
stantial amount of gas to forming galaxies, making them ca- 
pable of impacting the structural properties of galactic gas 
disks dSijacki et al.|2012| [Hobbs et al.|20T2 ). 



One way to improve the prospects of carrying out reli- 
able cosmological simulations without relying on substantial 
increases in available computational power is to improve the 
accuracy of the hydro solver for a fixed resolution or com- 
putational cost. For example, the primary reason for the 
required high resolution in Kaufmann et al. (2007) is to de- 
crease the importance of spurious hydrodynamical torques 
that occur at sharp density boundaries between dense galac- 
tic gas disk and the surrounding hot gas haloes. However, 
if one could remove the source of the spurious hydrody- 
namical torques and improve the thermal mixing properties, 
then it may be possible to relax the high resolution crite- 
ria to a more attainable level. It is possible that this could 
be achieved by either modifying the SPH algorithm (e.g. 
Ritchie & Thomas||2001| [ Pricel[2Q08l [Wadsley et al.|[2008l 
HeB k Springel 2010; Read e t al.|201Q[|Abel|2011| [Saitoh fc] 



disk-dominated Milky Way. However, Agertz et al. (2011) 
argued that - even though they were successful in repro- 



Makino^^2012^ ; Hopkins 2012j ) or by moving to a grid based 
code where these spurious hydrodynamical torques are not 
expected to occur (e.g. Okamoto et al.||2003 ). Indeed, in a 
recent study by Scannapieco et al. (2012) it has been shown 



for a single galaxy simulated at a high resolution via a zoom- 
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in technique that the choice for the adopted hydro solver can 
impact the galaxy morphology. 

In this paper, we explore the formation of gas disks 
in two cosmological simulations: one using a traditional 
density based SPH formulation as implemented in GAD- 
GET (Springel 2005) and one using a novel moving- 



mesh grid-based hydrodynamical solver as implemented in 
AREPO fSpringer'2010a"). Both GADGET and AREPO are 
massively parallel hydrodynamical simulation codes that use 
the same gravity solver and sub-grid physics, allowing us 
to isolate the impact of the hydro solver at an equivalent 
number of initial resolution elements and nearly equivalent 
computational cost. We study cosmological simulation runs 
with the two codes and find that the gas disk scale lengths 
associated with the cold gas disks formed in AREPO are 
systematically larger than their GADGET counterparts. We 
discuss the reasons for these differences including the im- 
pact of the numerical resolution. The work described here is 



an extension of the analysis presented in Keres et al. (2012) 



which discussed the properties of galaxies and haloes. In this 
work, we focus specifically on the properties of gas disks that 
form within the two simulations and compare their struc- 
tural properties. 

This paper is structured as follows. In Section 2, we 
summarize the numerical methodology. We then present our 
results for the gas disk properties in Section 3, and discuss 
the causes for the differences in Section 4. Finally, we sum- 
marize our findings in Section 5. 



2 METHODS 



We consider the simulations presented in |Vogelsberger et al.| 
(2011), which were performed with the cosmological codes 
GADGET ([Springel 2005) and AREPO (Springel 2010a). 
These simulations follow structure formation in an L = 
20Mpc box assuming a standard ACDM cosmology 

(Q^ = 0.73, Qo = 0.27, h = 0.7, as = 0.8). The simulations 
contain A^dm = 512^ dark matter particles and an equivalent 
number of gas resolution elements at the start of the simula- 
tion, giving a mass resolution of ttidm = 3.722 x 10^h~^MQ 
and TTigas = 7.444 x 10^h~^ Mq (note that the gas resolution 
element mass can change in AREPO due to mass advection 
across cell boundaries, but all cells maintain a mass within 
a factor of 2 of this target value). A comoving gravitational 
softening length of e = lh~^ kpc was used. 



2.1 Comparison Parameter Selection 

A central premise for our comparison is that GADGET and 
AREPO are very similar simulation codes which allows us 
to hold a large number of simulation parameters and proce- 
dures the same. In addition to containing the same grav- 
ity solver, the prescriptions for the radiative cooling of 



gas (Katz et al. 1996), the evolution of the ionizing back- 
ground radiation field ( [Faucher-Giguere et al. 2009), and 
star formation with associated feedback (Springel &: Hern-| 



quist 2003J are identical. The fundamental difference be- 
tween the codes lies in their handling of gas hydrodynamics. 
While GADGET uses a standard density based smoothed 
particle hydrodynamics (SPH) technique to evolve the gas, 
AREPO solves the fluid equations using an exact Riemann 



solver on a moving unstructured mesh based on a Voronoi 
tessellation. AREPO has a number of advantages over the 
standard SPH method due to its ability to, e.g., capture 
shocks more accurately and better resolve fluid instabili- 



ties ( |Springel||2010a||Sijacki et al.|2012] |Springel|2010b| ). In 

particular, AREPO can handle weak shocks more reliably 
than standard SPH codes which wash these features out 



because of the effects of artificial viscosity (see, e.g., [Keshet 
et aLpOO'sf . Also, with respect to Eulerian adaptive mesh re- 
finement (AMR) codes used in cosmology, AREPO employs 
a more accurate gravity solver (see, e.g.,'0'She a et al.|2005 ) 
and a spatial refinement that is continuous with the motion 
of the fluid. Although these traits have been demonstrated in 
isolated test problems ( Springel||2010a| Sijacki et al.||2012 ), 
we would like to understand how these numerical effects can 
impact the evolution of galaxies in cosmological simulations. 
Furthermore, because the same gravity solver and sub-grid 
physics prescriptions are used, we avoid some of the uncer- 
tainties that remain in code comparisons where these are 
allowed to vary (e.g., |Frenk et al.||1999| ). 

For our comparison, we choose to have both codes start 
from identical initial conditions and have the same initial 
number of resolution elements (A^dm = 512^, A^gas = 512^). 
From a practical standpoint, using the same number of ini- 
tial resolution elements leads to comparable computational 
cost for the two codes (with the AREPO simulation taking 
^ 30% more CPU time). This is an important considera- 
tion, because it is the CPU expense which sets limits on the 
size and complexity of simulations that can be run. However, 
since the two codes have roughly comparable CPU consump- 
tion for the same number of resolution elements, we simply 
note that neither code has a distinct advantage in this area. 

From a more physical standpoint, using the same num- 
ber of initial resolution elements results in similar mass res- 
olutions between the two codes. For the dark matter compo- 
nent where no inter-particle mass transfer is required, both 
simulations share identical mass resolutions. For the baryon 
component, the SPH particles in GADGET have a fixed mass 
in time while the cells in AREPO have a time dependent 
mass due to mass advection to their neighbors when solv- 
ing the Riemann problem across cell boundaries. We note, 
however, that we have included a refinement /de-refinement 
scheme that maintains all hydro cells in our AREPO sim- 
ulation within a factor of two of the SPH particle mass in 
our GADGET simulation ( Vogelsberger et al. 2011). As a 
result, both simulations share similar (although not iden- 
tical) mass resolutions. As discussed in 'Vogelsberger et al 



( |2011 ), we emphasize that the fixed particle mass in SPH 
is connected to inaccuracies made in solving the continuity 
equation, preventing this algorithm from correctly handling 
mixing. 

Since the dark matter components in both simulations 
evolve similarly (with some potential differences resulting 
from the influences of the baryon components), the same 
large scale structure and halo properties are present in both 



simulations (Vogelsberg er et al.||2011 ). As a result, we can 
compare gas disks which reside in a set of matching haloes 
that are identified as being present in both simulations. Fig- 
ures [l] and [2] illustrate this point for two matching haloes 
taken from our simulations. In these nested maps of the 
projected gas surface density, it can be seen that the dis- 
tribution of gas on large scales is similar in the two codes 
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Figure 1. Maps of the projected gas surface density for one object in the GADGET simulation at redshift z = 1. The central object 
has a halo mass of M = 2 x 10^'^h~^ Mq. Three nested views are shown to give a clear picture of the gas distribution over a large range 
of spatial scales. In the rightmost panel, the gas distribution around the central galaxy can be seen to be fairly clumpy and the galaxies 
themselves appear fairly compact. In this image, three galaxies are in the process of merging, which we have identified with red circles. 
It is helpful to directly compare this plot to Figure [2] which shows the same maps for the AREPO simulation. 
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Figure 2. Same as Figure^ but for the AREPO simulation. In contrast to Figure^ the gas distribution around the central galaxy in 
the rightmost panel has much fewer gas clumps. The three galaxies that are in the process of merging are highlighted with red circles. 



(see the leftmost panels of Figures [T] and [2]) . This is an ex- 
pected result, as the gas distribution traces the dark matter 
distribution on large scales. However, as one examines the 
distribution of gas on galactic scales (as shown in the cen- 
tral and rightmost panels), it becomes clear that there are 
prominent differences between the two codes. Specifically, 
while the gas in the GADGET simulation is distributed in 
a large number of compact and dense clumps, in AREPO 
gas has a much smoother distribution. Further inspection of 
Figures [l] and [2] shows that in the central region there are 
three galaxies - highlighted by red circles - in the process of 



merging. The effects of their mutual interaction can be seen 
in the case of AREPO, where tidal features are visible. Al- 
though these galaxies are present in GADGET as well, tidal 
features are much less prominent because the gas is more 
centrally concentrated and less rotationally supported. Per- 
haps the most striking finding from Figures [l] and |2] results 
from a comparison of the galaxy located nearest to the origin 
in these plots. While this object appears as a smooth, spa- 
tially extended disk in AREPO, the same object in GADGET 
is better described as a featureless blob. 

Figures [l] and |2] summarize the motivation for this work. 
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Figure 3. Projected gas surface density maps of 5 matched objects in GADGET and AREPO chosen at redshift z = with host halo 
masses ~ 10^2/^-1^0. There are clear diflFerences in the extent of the central gas disk. In addition, the prevalence of dense gas blobs is 
much higher in the GADGET simulation. 
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Figure 4. Projected gas surface density maps of 5 matched objects in GADGET and AREPO chosen at redshift z = 2 with host halo 
masses ~ 10^2/^-1^0. There are clear differences in the extent of the central gas disk. In addition, the prevalence of dense gas blobs is 
much higher in the GADGET simulation. 



Even though these simulations have been initiated from the 
same initial conditions, share the same feedback prescrip- 
tions, and use the same number of initial resolution ele- 
ments, the detailed morphological properties of the gas dis- 
tribution on galactic scales can be very different. This is a 
very important point, because it indicates that the hydro 
solver has a significant impact on the gas properties. While 
these differences can be fairly easily identified from the gas 
surface density maps, a primary goal of this paper is to pro- 
duce a detailed quantitative comparison of the sizes of the 
gas disks based on a large sample of galaxies matched be- 
tween the two simulations. 



2.2 Gas Disk Analysis 

In what follows, we contrast gas disks that form in the GAD- 
GET and AREPO cosmological simulations. To facilitate this 
comparison, we first identify a sample of matched galaxies 
from the two simulations. We start by building a catalog 
of all structure in the each simulation independently using 
the SUBFIND tool ( [Springel et aH|2001 |. We assemble a 
population of "matching" haloes by finding objects in the 
GADGET and AREPO simulations that have the same total 
mass to within 10% and potential minimum locations that 
are not offset by more than 25% of their half mass radii. 
We remove any pairs where the centre of mass is offset by 
more than 10/i~^kpc from the most tightly bound particle 
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Figure 5. Maps of the projected gas surface density for a typical matched galaxy in a Mnalo = halo. Red overplotted arrows 

denote the local gas velocity field. The AREPO galaxy (left panel) is significantly more rotationally supported than its GADGET 
counterpart (right panel). 



in either halo, as this may be an indication of a merging sys- 
tem. These selection criteria yield 1367 matching haloes in 
both simulations with total halo masses above 10^°/i"^ M© 
at redshift z = 0. Figures |3] and [4] show examples of the 
local gas density around a series of five such matching ob- 
jects in GADGET and AREPO at redshift z = and z = 2, 
respectively. 

In addition to showing differences in the spatial distri- 
bution of gas, the gas component of galaxies formed in GAD- 
GET and AREPO also differ kinematically. Figure [5] shows 
the maps of the gas surface density for an example matched 
galaxy with vectors indicating the local gas velocity field. 
Based on the topology of the velocity field, it is evident that 
the large gas disk in AREPO shows clear rotation about the 
disk's center. In contrast, the velocity field of the central 
galaxy for the same halo in the GADGET simulation, ex- 
hibits much lower level of circulation. 

To analyze the gas disk properties, we identify the grav- 
itationally bound cold and dense gas within each halo. We 
distinguish the diffuse hot halo from the colder, more dense 
rotating gas disk by making a cut in the T — p phase diagram 
at 



logi 



[K] 



6 + 0.251ogi 



lOio[M0 kpc- 



• (1) 



Our subsequent analysis depends on first removing the hot 
halo component before measuring the disk surface density 
profile. However, our results are not very sensitive to this 
particular cut in the phase diagram (i.e. moving the nor- 
malization of this cut up or down by a factor of 5 would not 
change our conclusions). The remaining cold and dense gas 
is translated to place the centre of mass at the origin and 
rotated to align the net gas angular momentum vector in 
the z direction. 

Before moving forward, we note that the definition of 



cold/dense gas cut defined in equation ([T]) can select mate- 
rial which is not part of the central gas disk. In particular, 
we find that high mass haloes in GADGET contain a popu- 
lation of low mass cold gas clumps which have no associated 
dark matter overdensity. The rightmost panel in Figure [l] 
shows an example of these gas clumps outside of the central 
object and a similar population of clumps can be found in all 
massive (i.e. IO^'^H'^Mq) haloes in the GADGET run. We 
note that such clumps are also seen in simulations of galaxy 
formation carried out with other SPH codes (e.g. Okamoto] 



et al.|[2008 Guedes et al. 2011). Most of these clumps are 
not part of the galactic disk - many of them have entered 
the halo for the first time, have not yet had any contact with 
the central galaxy, and are on very non-circular trajectories. 
However, the density and temperature of these clumps al- 
lows them to be selected as "disk mass" according to the 
definition outlined in equation ([T]). 

In principle, we could remove these cold gas clumps by 
imposing some additional criteria in our "disk gas" selec- 
tion, e.g., we could require disk gas to be on nearly circular 
trajectories or link together the central disk using a friends- 
of-friends (FOF) algorithm. Nonetheless, at this point we 
choose not to impose any additional selection criteria for 
two reasons. First, we find that de-selecting clump mate- 
rial can be sensitive to the details of the clump removal 
technique that we use. For example, while running a FOF 
algorithm on the cold and dense gas can efficiently remove 
clumps that are more than 15 or 20 kpc from the cen- 
tral galaxy, removing the most centrally located gas clumps 
can depend on our choice for the linking length. Second, 
since these clumps are only present in GADGET, we find 
that removing them tends to decrease the disk scale lengths 
obtained for massive GADGET galaxies without having any 
noticeable impact on the same objects in AREPO. Since one 
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Figure 6. Stacked mass profiles at z = for the low mass AREPO objects (top left), high mass AREPO objects (top right), low mass 
GADGET objects (bottom left) and high mass GADGET objects (bottom right). The true mass profiles are shown as solid lines, while 
the best fit exponential functions are shown as dashed lines. Each colour corresponds to an average surface density profile for all galaxies 
in a given halo mass range, with the average halo mass, Mavg, indicated in the legend. The vertical dot-dashed lines denote the spatial 
scale over which the gravitational potential is softened, i.e. 2.8 kpc, co-moving. 



of the conclusions in this paper is that gas disks formed in 
GADGET are indeed more centrally concentrated than those 
formed in AREPO, we have tried to avoid any steps in our 
analysis that could be perceived as artificially pushing us to- 
ward that result. Thus, in the following section on gas disk 
properties, we leave these clumps in our definition of disk 
mass and save additional discussion about their origin and 
impact for section [Xs] 



3 DISK COMPARISON 

To find the average surface density profile for gas disks, we 
stack the profiles for all objects in a given halo mass range. A 
representative set of stacked surface density profiles is shown 
in Figure |6] We find the best-fit exponential profile 



E(r) = So exp(-r/i?d) 



(2) 



for each surface density profile via a chi-squared minimiza- 
tion. Note from Figure |6] that most of the profiles can be fit 
using a single exponential. However, almost all profiles are 
distinctly steeper in the inner regions, and this is especially 
pronounced for GADGET galaxies. This compact central fea- 
ture is typically associated with a slowly rotating spheroidal 



component. Since we are primarily interested in the struc- 
ture of the spatially extended disks and because this inner 
region is not well-resolved in our current simulations (the 
softening length is lh~^ kpc co-moving, with gravity be- 
coming fully Newtonian after 2.8 kpc, co-moving, for the 
spline softening employed here (Hernquist & Katz 1989| )), 
we perform fits by both including and excluding the central 
2.8 kpc. We check the quality of each fit by comparing 
the integral of the best fit surface density profile in the fit- 
ted region to the true disk mass measured in the simulation. 
We note that when we exclude the central region all fits for 
both GADGET and AREPO return the cold, dense gas mass 
(as defined by equation |2| within 10%, indicating that our 
prescribed exponential functions are serving as good rep- 
resentations of the actual stacked surface density profiles. 
When we include the central region, the GADGET fits un- 
derestimate the disk mass by ~20-30% while the AREPO 
fits underestimate the disk mass by ^10%. This difference 
- which is more pronounced for GADGET - is caused by a 
concentration of material in the central region. 

The distribution of best fit parameters, So and i?d, of 
the exponential surface density profiles are shown in Fig- 
ure[7| We have performed the fits when including (left panel) 
and excluding (right panel) the central region. We find the 
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GADGET best-fit central densities are systematically low- 
ered when we exclude the central region, while the AREPO 
central densities changes less. This is a result of the fact 
that, on average, GADGET objects contain a larger fraction 
of their gas in the central region. More importantly, we find 
that the distribution of disk scale lengths is systematically 
different between the GADGET and AREPO galaxies, with 
the AREPO galaxies being described by larger scale lengths 
regardless of whether we account for the central region. This 
is a central result of this paper, and confirms the idea that 
the AREPO disks are larger based on visual inspection of the 
gas surface density maps. We emphasize that changing the 
number of galaxies in each stacked surface density profile, 
the exact normalization of our cold/hot phase boundary cut, 
or other detailed aspects of the analysis does not affect this 
conclusion. 

Figure [7| shows clearly the average offset toward larger 
disk scale lengths for the AREPO gas disks, we find it in- 
structive to show the best-fit exponential disk scale lengths 
as a function of host halo mass. Figure M shows the disk 
scale lengths (including the central 2.8 h~ kpc) as a func- 
tion of halo mass at redshifts z = 2 and z = in the left 
and right panels, respectively. This allows us to see more 
clearly how the gas disks from the two codes compare to 
one another at a fixed halo mass. For each bin we show disk 
scale length values obtained by first stacking the objects and 
then fitting the exponential surface density profiles (contin- 
uous lines), and by computing the median disk scale length 
from a set of individually fit galaxy surface density profiles 
(dashed lines). Hatched regions mark 25% and 75% of the 
distribution for the individually fit profiles. Regardless of the 
adopted procedure, the AREPO galaxies have larger scale 
lengths than their GADGET counterparts at both redshifts 
and for all halo masses. From Figure [S] it results that the 
AREPO disks are between 1.5 to 2 times larger than their 
GADGET counterparts. 

We note that for objects with a halo mass M > 
10^^'^h~^ Mq there is a discrepancy between the disk scale 
lengths obtained using our stacking procedure with respect 
to the median disk scale lengths obtained from individual 
objects in the GADGET simulations. The reason for this 
discrepancy is the presence of cold gas blobs surrounding 
the central gas disk. As discussed in the previous section, 
since these blobs are cold and dense they are included in 
our definition of "disk gas". Stacking many objects with a 
large population of blobs can then systematically bias the 
estimate of the disk scale length. Thus we note that the ap- 
parent trend of the disk scale length with the halo mass at 
the high mass end in the GADGET simulations as seen in 
the lower right panel of Figure [8] is affected by the presence 
of cold blobs. For these high mass systems, the central disks 
are just not well described by single exponential profiles. 
An analogous population of cold gas blobs is not present in 
the AREPO simulation. This is the primary reasons why the 
solid and dashed lines - representing the stacked best fit disk 
scale length and median of individually fit disk scale lengths 
- are in better mutual agreement for AREPO galaxies. 

Differences between the GADGET and AREPO galaxies 
can also be seen via a histogram of the gas specific angu- 
lar momentum, as shown in Figure [9] In this plot, where 
the specific angular momentum distribution for the example 
galaxies illustrated in Figure [5] is shown, there is a clear dis- 



tinction between the two codes. The AREPO galaxy shows a 
narrow distribution around jdrc with most of the gas on cir- 
cular orbits, and thus rotationally supported. The GADGET 
galaxy instead exhibits a much wider distribution around 
jz/ jcirc — 1 due to the significant gas component which 
is on highly non-circular orbits, with some material even 
counter-rotating, as evidenced by negative jz/ jdrc values. 

Figure [To] shows the specific angular momenta of the gas 
disks and of all baryons as a function of galaxy baryon mass. 
The hatched regions indicate where actual spiral and ellip- 
tical galaxies are located, as noted by Fall (1983). We have 



verified that the recent observational sample of |Courteau| 
et al. ( 2007 ) falls within the hatched region occupied by spi- 



rals. We find the specific angular momenta of the cold gas 
and stellar material in the AREPO simulation to be signifi- 
cantly larger than that of the GADGET objects. The larger 
specific angular momenta of the AREPO galaxies indicates 
that they are more rotationally supported, which accounts 
for their larger disk scale lengths. 



4 ORIGIN OF THE DISCREPANCIES 

In this paper we have presented a comparison of gas disks 



formed in cosmological simulations ( Vogelsberger et al. 



2011) performed with the SPH based code GADGET and 



the moving-mesh code AREPO. Both codes use an identical 
gravity solver and include the same physical processes (e.g., 
cooling, sub-grid model for star formation and feedback), 
but use fundamentally different hydro solvers. Whereas 
GADGET uses an SPH approach to evolve the gas, AREPO 
uses a finite volume scheme on a moving Voronoi mesh. 

Our primary conclusion is that the cold gas disks that 
form with AREPO are described by notably different sur- 
face density profiles than disks formed using GADGET. We 
showed that, on average, the cold gas disks in AREPO sim- 
ulations have significantly larger scale lengths compared to 
a matched sample of GADGET disks. Consistently, we also 
find higher specific angular momenta for the AREPO disks. 
Now that we have identified systematic differences in the 
disk scale lengths and angular momenta we discuss their 
principle numerical origins and address ways in which the 
discrepancies may be reduced. 



4.1 Spurious Hydrodynamical Torques 

It is well-known that conventional formulations of SPH suf- 
fer from artificial angular momentum transport at phase 
boundaries - like the hot halo cold-disk transition. For ex- 



ample, Okamoto et al. (2003) showed that SPH simulations 



are prone to angular momentum transfer at this interface 
and that this can have implications for disk formation in 



cosmological simulations. Moreover, Okamoto et al. (2003) 



demonstrated that shearing flows at phase boundaries are 
more accurately captured in grid-based hydro solvers. Thus, 
this is one particular area where we expect that the hydro 
solver included in AREPO should yield more accurate results 
compared to GADGET. 

One solution to this problem is to completely decouple 
the "hot" and "cold" particle neighbor searches (Okamoto 



et al.||2003||Marri k White 2003 ). In principle, this modifi- 
cation does a better job at estimating the local gas density 
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Figure 7. The best fit exponential surface density profile parameters to the whole disk (left) and with the central region excluded (right). 
The distribution of best fit parameters for and Eq is illustrated in the histograms. Clearly, AREPO objects have larger disk scale 
lengths than their GADGET counterparts. Including or excluding the central region does not alter this conclusion. 
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Figure 8. The best fit exponential disk scale lengths (in physical units) as a function of halo mass for GADGET and AREPO objects 
at 2; = 2 (right) and z = (left). Continuous lines indicate values for the stacked galaxies in each mass bin, while dashed lines denote 
the median i^^ obtained by fitting exponential profiles to individual objects. Hatched regions are 25% and 75% of the distribution. The 
AREPO disks are substantially larger at all halo masses and at both redshifts. 



using only particles that are in the same phase allowing for 
a cleaner separation of multi-phase gas boundaries. Since 
this will eliminate all hydro dynamical interactions between 
the hot and cold phase the spurious loss of angular momen- 
tum will be eliminated. However, decoupling the neighbor 
searches for particles based on their phase could also lead to 
artificial suppression of other physical phenomena that rely 
on direct interaction of multi phase gas such as ram pres- 
sure stripping or shock capturing (see, e.g., the discussion 



in |Marri Whitel|2003| for their procedure to address this 
issue) . 

Another solution proposed is to increase the simulation 



resolution substantially (Okamoto et al. 2003 Kaufmann 



[et al.||2007| . Increasing the resolution will reduce the influ- 
ence of the pressure gradient mis-estimation that occurs at 
density phase boundaries in standard SPH (see e.g. Agertz 



et al. 2007). [Kaufmann et al. (2007) tested the impact of 



resolution on idealized inside-out disk formation simulations 
using SPH and found that by increasing the number of SPH 



© 0000 RAS, MNRAS 000, 000-000 



10 Torrey et al 



particles in a halo above 10^ the artificial hydro dynami- 
cal angular momentum loss became subdominant to other 
torquing mechanisms. This problem is less severe (if present 
at all) in grid based codes like AREPO where both density 
phase boundaries and shearing flows are more accurately 



handled (Okamoto et al. 2003 Agertz et al. 2007 Sijacki 



et al. 2012). As a result, simulated gas disks in GADGET 



rapidly loose their angular momentum unless they have a 
very large number of resolution elements while the same gas 
disks are able to evolve without such severe angular momen- 
tum loss at the same resolution in AREPO. 

Since this artificial angular momentum transport is 
most prominent at low resolution, this is the main reason 
for the differences seen between the disk scale lengths of 
GADGET and AREPO objects at relatively low masses (i.e. 
< 10^^ Mq) and will contribute to the loss of angular 
momentum for higher mass systems (i.e. < 10^^h~^ Mq). 
We note that the most massive haloes considered in our cur- 



rent paper approach the resolution criteria set forth by |Kauf-"] 
mann et al.J (2007) (i.e. ^ 10^ particles per halo), so we 



do not expect spurious hydrodynamical angular momentum 
loss to be the dominant problem in our high mass simulated 
objects. However, for the low mass objects in our simula- 
tion which are by definition more poorly resolved, spurious 
angular momentum loss is bound to have a substantially 
larger impact. This conclusion is supported by examining 
the specific angular momentum content of the GADGET and 
AREPO galaxies - as shown in Figure ^] - which demon- 
strates that the discrepancy between the two codes is larger 
for low mass galaxies. 

The resolution dependence of spurious hydrodynamical 
angular momentum transport demonstrated in |Kaufmann| 
et al. ( |2007 t implies that increasing the resolution of our 
GADGET simulations by an appropriate factor (an increase 
in the number of SPH particles of 10"^ would give us the 
desired 10^ particles per halo at the low mass end) could 
improve the agreement between GADGET and AREPO. In 
fact, since massive galaxies are assembled via hierarchically 
merging smaller objects together, it is necessary to have > 
10^ SPH particles in all haloes - including low mass systems 
- to avoid spuriously loosing angular momentum in low mass 
systems that will ultimately become part of a well resolved, 
high mass galaxies (Kaufmann et al. 2007| ) . This is seemingly 
in accord with the results of very high resolution "zoom-in" 
simulations (e.g. Governato et al.|2004j that have found the 
angular momentum loss in gas disks formed in standard SPH 
simulations can be substantially reduced by increasing the 
particle resolution. However, while the 10^ SPH resolution 
elements per halo is a reasonable requirement for "zoom- 
in" simulations, this same requirement is not yet feasible for 
intermediate and low mass haloes in full cosmological box 
simulations. In that sense, we consider it an advantage that 
grid based codes such as AREPO do not suffer from this 
spurious hydrodynamical angular momentum loss even at 
resolutions well below 10^ particles per halo. 



4.2 Gas Heating and Cooling 

Although GADGET and AREPO share the same prescrip- 
tions for radiative gas cooling, non-adiabatic heat sources 
and mixing at phase boundaries can cause differences in the 
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Figure 9. The distribution of specific angular momenta in the 
disk for a typical matched galaxy in a Mnaio = 10^^ Mq halo. 
On the X-axis the ratio of the i-component of the specific angu- 
lar momentum to the expected specific angular momentum for a 
particle at that location on a circular orbit is shown. Both codes 
exhibit distributions that peak about 1, which corresponds to a 
rotationally supported gas disk. However, while the AREPO disk 
shows a narrow distribution, with most gas being on nearly cir- 
cular orbits, the GADGET object has a much wider distribution, 
with some gas on highly non-circular trajectories. 
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Figure 10. Specific angular momentum as a function of baryon 
mass for the gas (diamonds) and baryons (filled circles) for the 
matched sample of AREPO and GADGET galaxies. Hatched re- 
gions denote locations of spirals and ellipticals on this diagram 
as defined by pdl| flOSSt . 



growth of galactic gas disks. Vogelsberger et al. (2011) pre- 
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sented an analysis of the evolving thermodynamic gas prop- 
erties in our cosmological simulations. One conclusion from 
this analysis is that the cooling in AREPO is more efficient 
than in GADGET, which is driven primarily by differences 
in the cooling rates of diffuse halo gas. There are two main 
reasons for the differences in these cooling rates which are 
discussed in detail in Vogelsberger et al.| ( |2011 ) and summa- 
rized here. 

The first reason is differences in dissipative heating in 
haloes driven by the presence of turbulent energy. [Bauer <£] 
[Spr ingel' f2011') performed a comparison of the properties 
of simulated driven turbulence using GADGET and AREPO 
and found that while the two codes produce similar velocity 
and density power spectra for high mach numbers (i.e. for 
highly supersonic driven turbulence), there are prominent 
differences in the way turbulent power cascades to smaller 
scales in the subsonic regime. In AREPO, a Kolmogorov-like 
turbulent cascade is recovered which transports energy to 
smaller spatial scales. However, in GADGET, the turbulent 
large scale eddies are quickly dissipated close to the cooling 
radius and transformed into incoherent small-scale velocity 
noise which is converted into thermal energy as the velocity 
noise is damped out via artificial viscosity. This heats gas 
and inhibits cooling in GADGET haloes, driving part of the 
difference in the cooling rates seen in Vo gelsberger et al.| 
(20TT|. 



The second reason stems from differences in mixing be- 
tween the two codes, especially at density phase boundaries. 
For low mass galaxies, the gas cooling timescale in the halo 
is relatively short, such that the material is able to cool 
onto the central galaxy with similar efficiency in both codes. 
However, for a typical massive galaxy in our simulations, the 
diffuse halo gas becomes sufficiently hot and the gas cannot 
cool rapidly due to radiative losses. However, mixing that 
occurs around in-falling substructures as cool gas is hydro- 
dynamically stripped can substantially lower the local cool- 
ing timescale ( Marinacci et al.||2010 ). Idealized tests of gas 
stripping (A gertz et al.||2007| [Sijacki et aL]|2012[ ) show that 
this mixing will be suppressed in GADGET. 

A similar mixing boundary layer can develop at the in- 
terface of the central gas disk and the diffuse hot halo. This 
shearing phase boundary can generate mixing via Kelvin- 
Helmholtz instabilities which will be poorly resolved in GAD- 
GET dOkamoto et al.||2003t . Although the diffuse halo gas 
sitting just above the gas disk boundary may have long cool- 
ing times, continued efficient hot- mode gas accretion can be 
facilitated by mixing at this boundary. I Sijacki et al. (2012) 



demonstrated this point by examining the cooling rates of 
gas in an idealized gas sphere with GADGET and AREPO. 
They find that the cooling rates are nearly identical when 
no gas rotation is included, which validates that the cooling 
prescriptions in both codes are in fact functionally identical. 
However, when the same test is repeated with gas rotation 
the cooling rates for the two codes became discrepant, with 
AREPO cooling more efficiently. The more efficient cooling 
only sets in once a substantial amount of gas settles into a 
rotationally supported disk that can interact with the am- 



bient halo gas. ^Keres et al.| (|2012| present evidence for en- 
hanced cooling of this same sort in massive galaxies in the 
cosmological simulations by noting that the hot haloes in 
AREPO have cooler cores compared to GADGET and the 
gas in AREPO tends to be on radially inward trajectories - 



both of which are consistent with a scenario of less subsonic 
turbulence dissipation and of increased cooling in a mixing 
boundary layer between the diffuse halo gas and dense gas 
disk. 

The combined effect of "extra heating" in GADGET 
from the poorly resolved turbulent power cascade and "ex- 
tra cooling" in AREPO from better resolution of mixing at 
phase boundaries explain the global thermodynamic differ- 



ences seen m 



Vogelsberger et al. (2011). This also explains 



part of the differences in disk scale lengths that we see in 
this paper. Notably, high mass galaxies in AREPO will be 
more efficient at accreting material from the hot halo at late 
times which will help them maintain a gas rich disk. 

It is worth noting that several proposals have been put 
forward to improve the mixing at phase boundaries in SPH, 
which may bring the two codes into better agreement. For 
example, it has been suggested that by including a thermal 
conductivity term contact discontinuities and instabilities 
can be more accurately handled (e.g., Price|2008 Wadsley 
|et al .|l2008), which would improve the mixing picture with 
respect to the standard SPH. Alternatively, moving away 
from the density based formulation of SPH to an energy or 
pressure formulation can substantially reduce the artificial 
surface tension at contact discontinuities, which can help to 
resolve instabilities and mixing with higher accuracy (e.g.. 



Ritchie fc Thomas|2QQl| [Saitoh k Makino|2012| . In particu- 
lar, Hopkins (2012|) hasrecently shown how such a formalism 



can be derived from a variational principle, as in Springel] 
& Hernquist ( 2002 ) , resulting in a fully conservative version 



of SPH. Other changes, such as modifying the shape of the 
adopted smoothing kernel with a substantial increase in the 
number of neighbors used in the hydro calculations, or mod- 
ifying the momentum equation evaluated in the code, have 
all also shown promise in improving the performance of SPH 
at resolving instabilities (e.g., |Abel,2011, ,Read et al.|2010 ). 



4.3 Cold Clumps 

One important issue contributing to the differences in the 
disk scale lengths for high mass objects in GADGET and 
AREPO is the efficient accretion of low angular momentum 
gas via dense clumps. These clumps can be quite easily iden- 
tified using projected surface density maps of the gas around 
relatively massive galaxies residing in ^ 10^^ M© haloes in 
GADGET as is shown in Figure [3] and [4] at redshift z = and 
z = 2, respectively. In contrast, very few blobs are present 
in the AREPO galaxies aside from discrete objects which we 
have identified to be subhaloes with clearly associated dark 
matter components. 



4-3.1 Identifying Clump Material 

We can identify the clump particles by noting that they are 
overdense relative to the ambient hot halo density, colder 
than the ambient hot halo temperature, and do not have 
any substantial dark matter overdensity associated with 
them. To demonstrate this point. Figure [TT] shows a 2-D 
histogram of the material that is part of the primary sub- 
halo (as identified via SUBFIND) in the 5 matched haloes 
between lO^^/i"^ M© and 3 x lO^^/i"^ M© at redshift z = 
in density-radius space (also shown in Figure [3|. Although 
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Figure 11. The distribution of gas in the 5 matched haloes between lO^^ 3 >< iqIS/^-i M0 in AREPO (left) and GADGET (right). 
Red lines mark the boundaries we have chosen to separate the diffuse halo gas, central dense gas, and clump material are described in 
detail in the text. The fraction of gas that resides in the disk, /disk 5 halo, /halo 5 and clump, /dump regions are reported on each panel. 
Finally, the ratio of the clump gas mass to the disk gas mass, fc/fd is also shown. 



the clumps can be as or even more prominent at higher red- 
shifts, here we focus on this mass range and redshift in our 
analysis here for demonstrative purposes. 

We can place rough boundaries to break the density- 
radius space down into three physically meaningful compo- 
nents: central dense gas, diffuse hot halo gas, and cold/dense 
blobs. The central dense gas - which has been the subject of 
most of this paper - is concentrated into a relatively small 
region at small radii and has high density 



M 



< 0.1, P > lO^O l^vir 



• (3) 



The dense clumps are identified as dense material outside of 
the noted central dense gas region, i.e., 



clump 



= M 



> 0.1, p > 



10"^ 



(r/i?v 



Mq/i 
kpc^ 



(4) 



Finally, we assign all other material to be part of the diffuse 
hot halo which is characterized by relatively low density. 

The fraction of mass that resides in each region is noted 
on both panels of Figure [TT] In both GADGET and AREPO 
over 90% of the gas mass associated with these systems re- 



sides in the diffuse halo region. Vogelsberger et al. (2011) 



and [Keres et"ar ( |2012 ) analyzed the hot halo gas and found 



that more hot halo gas is present and that the temperature 
of this gas is hotter in GADGET compared to AREPO. This 
is consistent with our discussion from the previous section 
on the increased artificial heating in GADGET and increased 
mixing induced cooling in AREPO. 

The central dense region in AREPO contains about 8% 
of the total gas mass, compared to about 1% in GADGET. 
We identify two simple reason for this difference. First, more 
gas is able to cool into this region in AREPO for the rea- 
sons discussed in the previous section. Second, we find that 
the AREPO galaxies are able to maintain larger amounts 
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Figure 12. The distribution of clump masses is shown for GAD- 
GET (blue) and AREPO (red). We find that there is a substan- 
tially larger population of clumps in the GADGET simulation, 
and that these clumps build up around the mass resolution limit 
of our simulation which is marked by the vertical dashed line. 
We emphasize that the few "clumps" found in the AREPO sim- 
ulation are actually tidal features or recently stripped cold gas, 
which are identified as clump material according to our imposed 
density threshold. The overwhelming majority of the GADGET 
gas clumps do not contain dark matter and are not associated 
with infalling substructure or tidal features, as can be gleaned 
from Figures [1] and [3] 
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of gas in this region because the dense gas resides in a ro- 
tationally supported disk with intermediate star formation 
rates. In contrast, we find that the GADGET galaxies in 
these massive systems contain most of their gas mass in a 
very centrally concentrated region with high densities, which 
efficiently convert the gas supply into stars. 

We now turn to the "clump material" portion of this 
diagram. As we described earlier, material that resides in 
this region is far away from the central galaxy, but very 
dense compared to the ambient halo gas. Interestingly, we 
find a substantial population of cold clump gas in the GAD- 
GET objects that is not found in the AREPO systems. In 
terms of the fractional mass distribution, ^ 1% of the gas 
mass is in this region for the GADGET systems, while a 
negligible fraction is in the same region for AREPO objects. 
Furthermore, of the small gas fraction that does reside in 
this region for the AREPO objects, most of this material is 
associated with tidal features or recently stripped cold gas 
from infalling satellites. 

To quickly estimate the potential impact of these 
clumps on the growth of the central disk, we report the ratio 
of clump mass to disk mass in FigurepT] In GADGET we find 
that a comparable amount of mass is in the central dense re- 
gion as in the clump region. In other words, if these clumps 
are able to efficiently migrate toward the central object - 
which we will soon argue is the case - then they are capable 
of contributing a substantial amount of cold, low-angular 
momentum material to the central galaxy. We note that al- 
though these quoted fractional mass distributions depend on 
our choice for the boundary locations shown in Figure [TT] 
our conclusion that there is a substantial amount of mass 
in the clump region for the GADGET objects in this mass 
range is robust to any reasonable changes in the boundary 
definition. 



4.3.2 Clump Associations 

To classify these clumps, we select all overdense particles 
according to the cut presented in Figure [TT] and group these 
particles using a Friends-of-Friends (FOF) algorithm with 
a linking length of / = 3h~^ kpc yielding a set of particles 
belonging to each clump. Each clump then has a well defined 
position, velocity, and mass. 

One quantity of interest is the dark matter overdensity 
associated with these clumps. Since we have included only 
material that is bound to the primary halo within our FOF 
group, we expect that well defined substructure (e.g., in- 
falling dwarf galaxies) will have already been removed. We 
can verify that this is the case by checking the dark mat- 
ter overdensity associated with each clump, and comparing 
it to the spherically averaged dark matter density at the 
clump's position. In practice, we do this by finding the vol- 
ume associated with the N^^ nearest dark matter particle 
from the clump's center of mass, and compare that to the 
spherically averaged dark matter density measured in a thin 
shell at the same radius as of the clump. Although we have 
verified that this method would allow us to clearly identify 
substructure, we do not find any clumps in our sample with 
significant associated dark matter overdensities. We have re- 
peated our overdensity test using N*^ — 64, 100, 500 without 
any change in our results. 

The mass spectrum for these blobs is shown in Fig- 



ure [12] We find that the majority of the blobs contain just 
above or below 32 particles, which corresponds to the num- 
ber of nearest neighbors used in our SPH simulation as de- 
noted by the vertical dashed line in Figure [T2] Clumps tend 
to build up around the resolution limit of our simulation and 
as we have verified the mass spectrum changes accordingly if 
we increase or decrease the number of nearest neighbors. Al- 
though we do find a very small number of low mass cold and 
dense gas patches in the AREPO simulation, we note that 
these features are characteristically distinct from the cold 
gas clumps in GADGET. Specifically, these cold gas patches 
tend to be associated with tidal features or recently stripped 
gas from infalling satellites. So, while these cold gas patches 
meet the density based selection criteria that we have im- 
plemented and do not contain any clearly associated dark 
matter overdensity, we emphasize that their origin is very 
different from the large population of cold gas clumps that 
are seen the GADGET simulation. 



4.3.3 Clump Origin and Trajectory 

The origin of these clumps can be determined by tracing 
clump particles back in time. We find that the clump par- 
ticles originate in very mild gas over-densities within the 
filamentary structures of the intergalactic medium (IGM). 
As material in these filaments falls into the hot halo envi- 
ronment that surrounds central galaxies in massive objects, 
mild over-densities are amplified via hydrostatic pressure 
confinement. We find that the maximum past temperature 
of the clumps is well below the halo virial temperature, in- 
dicating that the blobs did not form via cooling instability 
of gas overdensities present in the hot halo (as studied in, 
e.g., Kaufmann et al.|2006 ), but rather remained cold during 
their accretion from the IGM. This point is demonstrated 
in Figure ^] which shows a 2-D histogram of the distri- 
bution of gas in radius-temperature space for 5 matched 



objects with halo masses just above M = 10 h~ Mq at 
redshift z — 0. As labeled within the plot the trajectory for 
clumps can be readily identified for the GADGET haloes, 
which allows cold material to migrate from large radii to the 
central object without ever heating substantially. No anal- 
ogous migration trajectory exists for the AREPO systems. 
The clump formation picture we have discussed here is con- 
sistent with the clump formation scenario outlined in [Keres] 
et al. (2009), where an analogous population of cold gas 



clumps were found originating from IGM filaments. Also in 
agreement with Keres et al. (2009), we find an accretion 



rate from these cold clumps in our GADGET simulation of 



O.5M0yr"' at redshift z = for M = lO'^/i"' 



systems. 

As a clump is gravitationally accelerated toward the 
halo center, the clump should begin to be disrupted and 
mixed via ram pressure stripping and the Kelvin-Helmholz 
and Ray leigh- Taylor instabilities. However, it is well known 
that these instabilities are poorly resolved and that ram 
pressure stripping is underestimated in the standard den- 
sity formulation of SPH. In particular, [Agertz et al.| ( |2007| 
and Sijacki et al. (2012) presented numerical experiments 
that showed cold blobs have artificially long survival times 
in SPH codes, while grid based codes like AREPO shred 
these clumps over substantially shorter timescales consistent 
with analytic expectations. As a result, mild over-densities 
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Figure 14. Maps of the gas surface density showing the trajectories of cold, dense gas clumps for a ~ l^^'^h~^ Mq halo in the 
GADGET simulation at redshifts z = {0.130,0.096,0.055} (left, central and right panel, respectively). A population of clumps (marked 
with diamond symbols) are identified in the leftmost panel and tracked forward in time. Dashed lines are denoting their trajectories. The 
vast majority of the clumps shown are moving toward the central galaxy on nearly radial trajectories, and are able to merge with the 
central galaxy. Thus, the clumps effectively provide a source of cold, low angular momentum material in the innermost regions, which is 
of entirely numerical origin. 



in the accreting filamentary material are condensed and frag- 
mented into a population of high density cold clumps that 
have substantially longer survival times in the GADGET sim- 
ulation compared to AREPO. Since idealized numerical ex- 
periments have shown that the survival timescale for cold 
gas blobs in GADGET is artificially long compared to ana- 
lytic expectations, we argue that the survival of these blobs 



in our GADGET cosmological simulations is a consequence 
of the same SPH deficiencies. 

When the clumps enter the virial radius for the first 
time, they have non-zero orbital angular momenta about the 
halo's center of mass that is consistent with other recently 
accreted material. However, as the clumps pass through the 
halo gas, they loose their angular momentum efficiently due 
to spuriously strong hydrodynamic drag forces ( |Tittley et al. 
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2001 ) and dynamical friction. This loss of angular momen- 



tum puts the clumps on increasingly radial trajectories and 
allows them to merge with the central object after only one 
or two orbits. Thus, the clumps become an efficient source of 
cold, low angular momentum gas feeding the central galaxy. 

Figure [M] specifically highlights the trajectories of a 
population of clumps in time. A set of cold, dense gas clumps 
are identified in the leftmost panel of Figure [M] and their 
subsequent trajectories are marked with dashed lines in the 
central and rightmost panels. By inspection of the marked 
clump trajectories, it can be seen that many clumps move 
on nearly radial trajectories and eventually merge with the 
central galaxy. Although each clump has a relatively low 



mass (- 10 - lO^'h'' M© - see Figure 12), a sufficiently 



large number of clumps fall into the central galaxy on a 
characteristic time scale of ^ IGyr. Thus, their cumulative 
mass amounts to a substantial fraction of the central galaxy 
gaseous mass. Furthermore, given that they arrive on nearly 
radial trajectories with low angular momentum, the clumps 
act as an efficient delivery source of low angular momentum 
fuel to the central galaxy. 

A similar population of blobs can be seen in the gas dis- 
tribution for relatively massive haloes (i.e. - lO^^h'^ Mq) 
in independent studies that used similar versions of GAD- 
GET (e.g., [van de Voort et al.|201l||Keres k Hernquist|2009 



van de Voort &: Schaye|20l'2 ). This is an interesting point be- 
cause the mass resolution used in van de Voort et al. (2011 ) 



is about a factor of 2 worse than the resolution used in our 
present study, and the resolution used in |Keres &: Hernquist| 
(2009 ) is a factor of 7 better. Both of these studies find that 



these clumps result from fragmented IGM filaments that are 
able to survive until they merge with the central gas disk. 
Despite their substantially higher resolution, Keres Hern-| 
quist (2009) find that these clumps tend to form just above 
the resolution limit of their simulation - implying a signifi- 
cant change in the clump mass spectrum from what we have 
found here. 



Recently, Hobbs et al. (2012) have made use of a new 
flavor of SPH (Read k Hayfield "2012) to study the forma- 
tion and impact of "blobs" in standard SPH simulations and 
found their origin to lie in artificial thermal instabilities that 
can originate from a small number of particles with low en- 
tropies with respect to their neighbors. Since the standard 
formulation of SPH in GADGET lacks any inter-particle fluid 
mixing, these artificial thermal instabilities are allowed to 
grow and ultimately form a population of dense gas clumps. 



Hobbs et al. (2012) show that this artificial thermal instabil- 



ity can be averted by including thermal conductivity (e.g., 
[Pri ce 2008). This issue does not arise with AREPO because 
cells exchange entropy with their neighbors when mass is 
advected across cell boundaries, which results in a physi- 
cally motivated homogenization of the fluid entropy. The 
presence of these clumps in independent numerical studies 
of varying resoluti on (e.g. [van de Voort et al.||2011| [Keres[ 
\k Hernquist^^2009 van de Voort k Schaye 2012| seems to 
imply that these clumps will continue to form and survive at 
higher resolution even though the specific properties, such 
as the mass spectrum, will change unless modifications are 
made to the standard SPH hydro solver. 



5 CONCLUSIONS 

Here we presented a comparison project aimed at study- 
ing the properties of gas disks formed in cosmological sim- 
ulations performed with two very different hydrodynami- 
cal codes - GADGET and AREPO. Our comparison started 
from an identical set of initial conditions which were evolved 
forward in time with the two codes which have the same 
gravity solver and while holding fixed the initial number of 
resolution elements, and radiative cooling, star formation 
and feedback prescriptions. However, GADGET and AREPO 
adopt very different approaches for solving the hydro dy- 
namic equations. While GADGET uses a standard density 
formulation of SPH, AREPO solves the fluid equations using 
a Riemann solver on an unstructured moving mesh. In many 
ways, the hydro solver included in AREPO has accuracy ad- 
vantages over the SPH solver used in GADGET which can 
be clearly demonstrated in, e.g., idealized shock tube tests, 
driven turbulence tests, and hydro instability tests. The goal 
of our comparison was to understand the impact of the hy- 
dro solver on the formation of gas disks in fully cosmological 
simulations at a comparable resolution. Our primary conclu- 
sions are as follows: 

• After fitting the gas disks with best-fit exponential sur- 
face density profiles, we find that the AREPO gas disks are 
systematically larger than their GADGET counterparts. This 
corresponds to gas disks in GADGET having lower specific 
angular momentum compared to the matching set of disks 
formed in AREPO simulation. 

• The primary reason responsible for the differences in 
gas disk scale lengths between the two codes changes as a 
function of the number of resolution elements and physical 
environment of the host halo. 

• For low mass objects, low resolution leads to spurious 
angular momentum transport from the cold disk to the dif- 
fuse hot halo in the GADGET simulation. This spurious an- 
gular momentum loss is a well-known and documented issue, 
which can be alleviated by moving to increasingly higher 
resolution in test problems or "zoom-in" simulations. How- 
ever, for large cosmological box simulations, like those we 
have presented here, the resolution needed to suppress this 
spurious angular momentum loss is not yet attainable. Grid 
based codes - such as AREPO - are not expected to suffer 
from this same problem and can therefore provide a more 
accurate answer for the same number of resolution elements 
and comparable CPU time. 

• Poorly resolved subsonic turbulence in GADGET results 
in dissipative heating of the gas near the cooling radius. This 
inhibits the accretion of gas onto the central galaxy. In the 
presence of turbulent energy, AREPO correctly recovers a 
cascading Kolmogorov-like power spectrum, so no analogous 
artificial heating source is present. 

• For high mass objects, the cooling timescale of hot halo 
gas can become comparable to the Hubble time which ef- 
fectively shuts off fresh gas accretion. However, mixing at 
density phase boundaries - such as the interface between 
the cold gas disk and the hot halo - can substantially in- 
crease the gas cooling rates. This allows for more gas to 
cool onto the disk in AREPO given that the mixing at phase 
boundaries is suppressed in GADGET. 

• For high mass objects, the efficient delivery of low an- 
gular momentum gas in the form of cold gas clumps causes 
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the central gas disks in GADGET to be much more centrally 
concentrated than in AREPO. These clumps form from frag- 
mented IGM filaments and rapidly migrate to the potential 
minimum as they loose their angular momentum to hydro- 
dynamic drag against the ambient hot halo. The absence of 
these blobs in AREPO is attributed to the efficient disrup- 
tion of clumps via ram pressure stripping and the Kelvin- 
Helmholtz and Ray leigh- Taylor instabilities - all of which 
are poorly handled in GADGET. 
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